In [1]:
%matplotlib inline
from common import *
datadir = os.path.join("//media", "disk", "Data")
#datadir = os.path.join("..", "..", "..", "..", "..", "Data")
In [2]:
south_side, points = load_data(datadir)
In [3]:
masked_grid = grid_for_south_side()
print("Total number of cells in region:", masked_grid.xextent * masked_grid.yextent)
print("Number of cells which intersect geometry:", np.sum(~masked_grid.mask))
In [4]:
fig, ax = plt.subplots(ncols=2, figsize=(16,8))
for a in ax:
a.add_patch(descartes.PolygonPatch(south_side, fc="none", ec="Black"))
a.add_patch(descartes.PolygonPatch(south_side, fc="Blue", ec="none", alpha=0.2))
xmin, ymin, xmax, ymax = south_side.bounds
a.set(xlim=[xmin-500,xmax+500], ylim=[ymin-500,ymax+500])
pc = open_cp.plot.patches_from_grid(grid_for_south_side())
ax[0].add_collection(matplotlib.collections.PatchCollection(pc, facecolor="None", edgecolor="black"))
pc = open_cp.plot.patches_from_grid(grid_for_south_side(xoffset=120, yoffset=120))
ax[1].add_collection(matplotlib.collections.PatchCollection(pc, facecolor="None", edgecolor="black"))
ax[0].set_title("Original gridding")
ax[1].set_title("Offset by 120 meters north and east")
None
In [5]:
cts_pred = retro.RetroHotSpot()
cts_pred.data = points
cts_pred.weight = retro.Quartic(bandwidth = 1000)
cts_risk = cts_pred.predict(end_time = np.datetime64("2011-09-28"))
cts_grid_risk = open_cp.predictors.grid_prediction(cts_risk, masked_grid)
grid_pred = retro.RetroHotSpotGrid(grid=masked_grid)
grid_pred.data = points
grid_pred.weight = retro.Quartic(bandwidth = 1000)
grid_risk = grid_pred.predict(end_time = np.datetime64("2011-09-28"))
grid_risk.mask_with(masked_grid)
In [6]:
fig, ax = plt.subplots(ncols=2, figsize=(16,8))
for a in ax:
a.set_aspect(1)
a.add_patch(descartes.PolygonPatch(south_side, fc="none", ec="Black"))
ax[0].pcolormesh(*cts_grid_risk.mesh_data(), cts_grid_risk.intensity_matrix, cmap=yellow_to_red)
ax[0].set_title("Continuous prediction")
ax[1].pcolormesh(*grid_risk.mesh_data(), grid_risk.intensity_matrix, cmap=yellow_to_red)
ax[1].set_title("Grid prediction")
None
In [8]:
evaluator = open_cp.evaluation.HitRateEvaluator(RetroHotSpotEval(masked_grid, points))
evaluator.data = points
result = evaluator.run(time_range(), range(0,51))
In [9]:
frame = to_dataframe(result.rates)
fig, ax = plt.subplots(figsize=(16,6))
(frame[5] * 100).plot(ax=ax)
(frame[10] * 100).plot(ax=ax)
(frame[15] * 100).plot(ax=ax)
ax.set_ylabel("% hit rate")
ax.legend(["{}% Coverage".format(x) for x in [5,10,15]])
None
In [10]:
frame.describe()
Out[10]:
In [11]:
fig, ax = plt.subplots(figsize=(6,6))
d = frame.describe()
ax.plot(d.columns, d.ix["mean"]*100)
ax.plot(d.columns, d.ix["25%"]*100)
ax.plot(d.columns, d.ix["50%"]*100)
ax.plot(d.columns, d.ix["75%"]*100)
ax.legend()
ax.set(ylabel="% hit rate", xlabel="% coverage")
ax.set_title("Hit rate against coverage")
ax.set(xlim=[0,50], ylim=[0,100])
None
In [12]:
predictions = []
times = list(result.details)
times.sort()
for time in times:
predictions.append( result.details[time].prediction )
In [ ]:
from matplotlib import animation
matplotlib.rc('animation', html='html5')
In [ ]:
fig, ax = plt.subplots(figsize=(8,8))
ax.set_aspect(1)
ax.add_patch(descartes.PolygonPatch(south_side, fc="none", ec="Black"))
In [ ]:
def animator(n):
pred = predictions[n]
ax.collections = []
mesh = ax.pcolormesh(*pred.mesh_data(), pred.intensity_matrix, cmap=yellow_to_red)
ax.set_title("Risk density on {}".format(times[n]))
return mesh
In [ ]:
animation.FuncAnimation(fig, animator, frames=len(predictions), interval=100, blit=False)
In [7]:
# This takes a long time to run..
with open_cp.pool.PoolExecutor() as executor:
futures = []
for xoffset in range(0, 250, 5):
for yoffset in range(0, 250, 5):
masked_grid = grid_for_south_side(xoffset, yoffset)
task = RHS_Eval_Task(masked_grid, points, (xoffset, yoffset))
futures.append( executor.submit(task) )
grid_effect = {}
for key, result in open_cp.pool.yield_task_results(futures):
grid_effect[key] = result.rates
In [16]:
import pickle, lzma
with lzma.open(os.path.join("grid.pic.xz"), "wb") as f:
pickle.dump(grid_effect, f)
In [20]:
import pickle, lzma
with lzma.open(os.path.join("grid.pic.xz"), "rb") as f:
grid_effect = pickle.load(f)
In [21]:
def to_mean(result):
frame = to_dataframe(result)
return frame.describe()[5]["mean"]
means = { key : to_mean(value) for key, value in grid_effect.items() }
In [23]:
first_axis_keys = set(k[0] for k in means.keys())
second_axis_keys = set(k[1] for k in means.keys())
frame = pd.DataFrame({k:{l:means[(k,l)] for l in second_axis_keys} for k in first_axis_keys})
frame
Out[23]:
In [30]:
data = np.asarray([np.asarray(frame[x]) for x in frame.columns])
data.shape
Out[30]:
In [43]:
fig, axes = plt.subplots(ncols=2, figsize=(17,6))
ax = axes[0]
ax.hist(data.flatten(), bins=np.linspace(0.195, 0.23, 20))
ax.set_title("Histogram of mean hit rate")
ax.set_xlabel("Man hit rate")
ax.set_ylabel("Count")
ax = axes[1]
mappable = ax.pcolor(data)
ax.set_title("Mean hit rate for offset of grid")
ax.set_xlabel("X offset, times 5m")
ax.set_ylabel("Y offset, times 5m")
cbar = plt.colorbar(mappable)
cbar.set_label("Mean hit rate")
None
In [ ]: